Equation of State and Phases of Polarized Unitary Fermi Gas 
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The equation of state of the partiaUy polarized two component Fermi gas at zero temperature 
in the unitary limit is computed by ab initio auxiliary field Monte Carlo method. From this, we 
obtain the critical ratio of the chemical potentials at the phase transitions. The value of 

fill 11^ at the transition between the fully paired superfluid and the partially polarized phases is 0.11 
while the critical value at the phase transition between the partially polarized phase and the fully 
polarized normal fluid is —0.59. We also determine the radial boundaries of the phase transitions of 
the Fermi gas in the harmonic trap as function of the total polarization. We find that beyond the 
critical polarization 0.65, the fully paired superfiuid core disappears in the trapped Fermi gas. 

PACS numbers: 03.75.Ss, 05.30.Fk, 02.70.Uu, 21.65.-f 



Dilute fcrmion gases such as those of ^Li and "^''K 
are quantum mechanical systems w\t\v controllable short 
range and strong interactions. They offer an ideal 
test bed for our knowledge of the quantum many-body 
physics. The properties of these Fermi gases can be 
probed experimentally The interaction is described 
by the dimensionless parameter askp where Os is the s- 
wave scattering length and A:^^ is the Fermi momentum. 
In the weakly interacting (so-called BCS) regime where 
l/uskp << and up to the strongly interacting regime 
with l/askp sa (also called unitary regime), one or 
more non-trivial phases have been suggested [1, 0, [E| 
for the Fermi gases with spin imbalance. These polar- 
ized atomic gases hold resemblance to the case of magne- 
tized superconductivity [6,]. Here, instead of the external 
magnetic field, we assume unequal chemical potentials. 
The constraints given on the chemical potentials of the 
different fermion species suggest the existence of one or 
more intermediate polarized phases However, these 
phases are hard to study theoretically since the mean 
field approaches are not quantitatively accurate, and the 
numerical techniques such as the Fixed Node Diffusion 
Monte Carlo (FN-DMC) method requires the knowledge 
of the physically motivated guiding functions in the first 
quantized form. 

The system we consider in this article is that of the ide- 
alized Fermi gas consisting of two^ ,|-spin) species with 
equal masses. We assume control on each one of the 
chemical potentials (/i|,/ij^) and the physically measur- 
able quantities are the densities n| and nj^ (0 < < ^t)- 
In the spin symmetric phase, the existence of the gap 
is manifest in the fact that the densities are not sen- 
sitive to a small difference of the chemical potentials 
S/i = (/Z| — /ix)/2. The superfluid phase imposes the 
constraint 6fj, < A [8]. Here A is the usual superfluid 
pairing gap of the symmetric system. This condition 
gives the lower bound on the critical y = Hi/fJ-^ deflned 
as Yi 0. We use a capitalized notation to indicate 
the upper or lower bounds while the lower case notation 
Ux corresponds to the actual critical value at the phase 
transition. At a specific value of yi > Yi (or correspond- 
ingly at a critical (5/i), the fully paired superfluid (SF) 



undergoes phase transition into the partially polarized 
phase {PP, see Fig [3]) and the densities become unequal 
(n| < rij). Recently, two possibilities were considered. 
In the first, the SF phase could transition into the po- 
larized normal phase going through a phase separated 
mixture of the superfiuid and the partially polarized nor- 
mal fiuid 0, . This transition that is assumed to be of 
the first order, is characteristic of the weakly interact- 
ing regime and also of the unitary limit. Another possi- 
bility is that the SF phase undergoes the second order 
phase transition and becomes a homogeneous polarized 
superfluid that accommodates the excess of one species. 
This was suggested for small polarizations in the unitary 
limit by Carlson et al. [l3|. This phase is alternatively 
called gapless or polarized superfluid phase {SFp). The 
recent work by Pilati et al. [11[ considers both possibil- 
ities and suggests that the gapless homogeneous phase 
may occur in the l/uskp > regime for moderate po- 
larization. In the limit of the complete polarization, the 
system is in the normal fully polarized(iVi?p) phase with 



Ni spin up particles and /i| 



> 0. This 



system is also insensitive to changes of Sfi as long as 
fi^ << 0. When ni > the energy difference between 
7V| -f- 1| and systems, the system phase transitions 
into the partially polarized(PP) phase. This defines an 
upper bound forw known as lo Several authors have 
shown 0, H, [i2j [3, that a simple variational solution 
of non-interacting Npp + interacting impurity gives an 
upper bound Yq reasonably close to the actual threshold 
value yQ. 



In this letter, we construct the equation of state of 
the unitary Fermi gas connecting the limits x = — — 

and X = 1. Then we estimate the actual yi (> Yi) and 
2/0 (~ Yq) as a direct application of the knowledge of the 
equation of state. We also verify the consistency with the 
previously reported values of Yo and Yi . We also present 
the equation of state in terms of the grand canonical po- 
tential (pressure) and the density profiles of the trapped 
gas. In order to do so, we implement the canonical en- 
semble auxiliary field Monte Carlo (AFMC) formalism 
at zero temperature. AFMC is usually formulated in the 
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second quantized form. In principle, it does not depend 
on the particular choice of the basis. It can be applied 
to the finite temperature , [13, as well as zero tem- 
perature 113 Fermi systems. With the cost of in- 
troducing a set of additional integration variables, the 
time propagator can be expressed in the basis set of one 
particle orbitals. Then, the multi-dimensional integra- 
tions over the additional auxiliary variables are carried 
out by the Monte Carlo method. The decomposition of 
the attractively interacting potential into negative eigen- 
values avoids the essential sign problem associated with 
the complex time evolution matrix. However, for the 
imbalanced systems the sign problem appears for long 
enough time evolution for which we give a simple prac- 
tical solution. We briefly outline the zero temperature 
canonical formalism with fixed particle numbers and 
iVj. We assume a zero range interaction of strength g 
between the particles of different spin species. We have 
no exchange interaction and the Hamiltonian adopts the 
form 



(r) 



(1) 



where ^<T(r) and ^t(r) are the usual fermion field oper- 
ators. We implement the solution of this Hamiltonian in 
a cubic volume of lattice sites. In this case, the con- 
stant g is lattice renormalized coupling by the cutoff kc in 
the momentum space: 1/g = 1/gb — l/^^X]|k|=o ^/(^Ck)- 
Here, the bare coupling constant gb = Anh^as/m, ek = 
?l^fc^/(2m), and fl is the volume of the system. 

The time evolution operator can be decomposed into 



-Ar-HlM 



a product of smaller steps e" 

Jg-Arr/2g-ArVg-ArT/2jM ^^^^^ ^ ^ ^ rj.^^ 

tributions of the one-body operators such as e^'^'^^/^ 
(with T = — X](T can be easily estimated 

in the one-body basis (that is, in the momentum space 
connected by Fourier transform). However, the interac- 
tion term (e"'^'^^) needs to be treated before it becomes 
computable. Since the operator ricr(r) = ^'3^(r)*I'o-(r) has 
eigenvalues or 1, n^(r) = ^^(r) and we can write [lB| 

"tW^iW = i(nT(r)+"i(r))'-^(riT(i-) + ^i(r)) • (2) 

The last parenthesis is a sum of one-body operators. 
However, the first parenthesis is a square of one-body 
operators. We introduce a set of one dimensional contin- 
uous variables at each lattice site r and use the Hubbard- 
Stratonovich transformation to get 

-T^E9("T(r)+fti(r))" 
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FIG. 1: (Color online) Imaginary time evolution for different 
lattice volumes Nj^ and particle numbers A'' = 2A'^| = 2Ni^. 
We assume h = 1. In general, the convergence is reached 
after tEfg > 3. We can see that for iV; > 7 and iV > 38, the 
ground state energy per particle converges to ~ 0.44 in units 
of Efg- 



Then the time evolution operator can be expressed as 
- ' V[x]G[x]U[x] (4) 



measure T^[x] — 
< M) includes the aux- 
in the time slices {Tm}- 
auxiliary 



the 



variables: 



where the integration 
nr,mC?3:(r,r„)^ (0 < m 
iliary variables {x(r,T,„)} 
G[x] is the Gaussian of 

G[x] — e Discrete auxiliary fields 

can also be used as shown in the references [igI . [isj . The 
Gaussian factor G[x] can be sampled directly leaving 
only the integrand U[x] that is used as the probability 
density for the Metropolis random walk. From the 
stability argument we know that Eq > 0, thus we have 
the bounds d{n{T))/dT < and d^{n{T))/d^T > 0. 
The r is pushed to the limit tq where the plateau with 
(9(7i(r))/9r < is reached (Fig[T|). Then we enter into 
the sign problem region where further evolution results 
in alternating signs of the spin unmatched {Nij ^ N^) 
fermion propagator. Instead of taking samples at 
various points of r > tq, we perform separate runs with 
evolution up to tq to get the samples. We determine 
that the convergence to the ground state occurs after 
tqEfg > 3 (?i set to 1). 

In the zero temperature formalism, we assume an ini- 
tial wave function that is not orthogonal to the ground 
state The ground state is projected out by taking 
lim e^'^^y&f. This is analogous to the DMC. However, 

T — >00 

here *t is a Slater determinant represented by iVj, x iV 
matrix where Nt, is the size of the basis set and N the 
number of particles. For the general quasiparticle creator 



J2j with a\ = plane wave creation operator. 



{Dji} are the elements of the matrix 



representing 

the state ni=i n In our case, the initial state is 

constructed by completely filling the lowest (and Ni ) 
plane wave states with equal amplitude {Dji = Sji). The 
successive applications of the short time evolution oper- 
ator yield a matrix of the same form and dimension. We 
can separately treat the different spin components of the 
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FIG. 2: (Color online) Equation of state of the unitary Fermi 
gas as function of the relative density x. We define the energy 

Efg^ = 'i^-^^^rl^r—- The AFMC results were obtained for 
A'^l = 19 in a 7^ volume (circle) . For the comparison purpose, 
we also plot the FN DMC resuhs of PP (triangle up, Ref [§]) 
and SFp (square, Ref [lC|).The fit to the equation of state is 
consistent with the convexity of the thermodynamic potential. 



wave function, SI/q = C/| (r)*^,! C/| (r)\I>t^l and have 
the density operator in the momentum space 



it 



(5) 



Here Uair) represents the spin a time evolution oper- 
ation and *o,(T = Ua-{T)'ift,a- The energy can be cal- 
culated as (n) = E f^k(ak,^ak,<T> + .gI](f^T(I■)"-i (!■))• 



For ^ Ni and large 



The sampling probability is 

enough r this probability can be negative. Thus, we 
sample instead |(4'f|^o)| and the normalization becomes 
samples I *o) / 1 (*t I *o) I • This sigu problem may arise 
because *t,T and ^t,i can be different when 7^ -^j.- 
We notice that for the case considered here where = 
19 and tEfg ~ 3, the sign overlap of the determinants 
is still positively biased and the normalization non-zero. 
This is true even in the extreme polarization of -I- 1| 
case. 

The possibility of the first order phase transition be- 
tween SF and Npp along < a: < 1 without any in- 
termediatepartially polarized phase(PP) was discussed 
by Cohen [8|. However, from the constraints for Yq and 
Yi that establishes a strict inequality Yq < Yi this 
is not a likely scenario. Hence, we assume the exis- 
tence of PP phase and the first order phase transition 
SF ^ PP. Then, we try to identify qualitatively dif- 
ferent regions of E{x). We identify that Xc ~ 0.42 sep- 
arates those regions (see Fig [5]). At x ^ 1, apparently 
there is a hump in the energy which we interpret as the 
consequence of the finite size of the system. It would be 
mostly due to the finite size pairing gap that will scale 
as ~ A/N^. In the range of a:c < < 1, we fit the 
equation of state (data points at a; = Xc and 1) with 
the functional form f{x) consistent with the convexity 



constraint and the Maxwell construction 7, 's'l . Thus we 
take the form f{x) ^ {a + bxf/^ (see EqlH). All the 
data points in < x < were best fitted by a polyno- 
mial function of third power. At x = 1, E / {N Epc) = 
0.44(2) (with EpG = | ^'(3^'(^t+'h))'^^ )j^ close match 
with the previously known values 13, SH. This num- 
ber is free of any sign errors and convergent at different 
system sizes and filling factors we have considered (see 
Fig [J). The one particle chemical potential measured as 
{E{N^, 1|) - E{N^,Q{))/EpG] gives -0.99(2)consistent 
with the variational calculations 0, 0, [13, [3, [l3| and 
somewhat higher than Prokof'ev and others' diagram- 
matic Monte Carlo estimate of -1.03 [l5|. Our value 
corresponds to ~ = —0.59. On the other hand, 
the value of yi calculated by taking at x ^ 1 is 

yi « 0.11, well above the lower bound Yi « —0.1 and 
Y. =j2Cf^^ — 1 w —0.1. Yc was defined in the references 

[l3| as the critical value for y where if Yq = Y^ = Yi 
the PP phase would disappear. Clearly, this is not the 
case according to our results. Experiments seem to give 
a rather wide range of possible values for the quantity 
7 = (1 - yi)/(l - yo) = (1 - rl)/{l - rg (see Fig H for 
the definitions of rn and ri and the Ref Q for that of 7). 
7 = 0.70 in Ref ^ and 7 = 0.56 in Ref [H]. In these 
experiments the effects of the finite temperature and the 
expansion introduce additional corrections. From our re- 
sults of j/o and yi we get 7 « 0.56. At the SF ^ PP 
phase transition Sfi/Epc — |f(l — yi)/(l -I- j/i) ~ 0.58 
where ^ is the energy per particle in the SF phase in 
units of EpQ. Thus, the superfluid pairs start to break 
at S^/A > 0.70 (using A/Epc = 0.84(4) from reference 
[l0| ) while the completely polarized phase is reached for 
(5/i/A > 3.37 (see Fig[3l). This picture is consistent with 
our earlier assumption that along the x direction regions 
of pure PP and mixture PP + SF exist rather than 
Npp -\- SF mixture in the whole partially polarized re- 
gion < X < 1. 

The free energy density of the polarized Fermi gas can 
be written as a function of the partial densities £{nij, n^), 
where remains fixed and changes from to rij . The 
chemical potentials are given by = dE/dua- Alter- 
natively, the pressure can be obtained by the Legendre 
transform V{^J,^,fJ,i) = /i|ri| -I- /ij^nj^ — £{n^,ni) where 
Ua- = dP/dfia- From the dimensional argument, we can 
write the energy and the pressure in terms of the dimen- 
sionless functions /(x) and h{x) fi^,^ 
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[^^My)?^' ■ (6) 



Directly from the Monte Carlo output, we can give the 
equation of state either in the £ vs (see Fig [2]) or V vs 
Ha {see Fig [3]) format. The h{y) vs y figure (Fig [3]) shows 
the pure phases in terms ofy. In this figure, we can easily 
locate the relevant values of y. The inset figure for the 
pressure shows discontinuity in the slope at the point of 
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FIG. 3: (Color online) We can clearly locate the values of 
yo r; Yq, Yi and yi{> Yi) in this figure for h{y) (Eq[6l. 
The discontinuity in h'{y) (correspondingly in the densities) 
at yi is due to the first order nature of the phase transition 
between SF and PP. For y < yo, we have Npp phase, for 
yo < y < yi the PP phase and for y > yi fully paired SF 
phase. Here we assume yo ~ Yq. The inset figure corresponds 
to the pressure in units of the superfluid pressure Psf as a 
function of Sfi/ A (we take A/Efg ~ 0.84 [l3|). Here we set 
(/Xf + Hi) 12 constant. The plateau in the pressure for small 
values of Sji/A is due to the existence of the pairing gap. It 
is the projection of the straight segment of h{y) at y > j/i. 




FIG. 4: (Color online) Local polarization a{r) at different 
total polarizations Pol. The radius r is given in units of r„ac. 
Tvac is defined as A| = y(r„ac). Pol is controlled by the ratio 
of the global chemical potentials A|/A| . At the critical Pole ~ 
0.65 we have the threshold where ri > appears (at smaller 
values of Pol), ri is the boundary of the fully paired core 
with cr(r]~) — and o-{rf) > 0. At smaller polarization(Poi), 
we can see the jump in the local density at ri corresponding 
to the first order phase transition SF <-> PP. The radius ro 
equals the position where a{ro) = 1. This transition is of the 
second order and the a{r) behaves continuously. 



SF PP phase transition while the slope is continuous 
at the transition PP ^ Npp. 



For the case of the trapped gas, the local density ap- 
proximation (LDA) defines the local chemical potentials 
/io.(r) = Act — V{r) with A^ = global chemical poten- 
tial and V{r) — imw^r^, the harmonic trapping po- 
tential. Using the equation of state and the thermody- 
namic relations, we can establish the local density for 
each spin species Uair) = ncr(/i| (r), /i| (r)) in terms of 
h{y) and h'{y). We can define the local polarization 
a(r) = (rij (r) — n|(r))/(n| (r) -I- (r)) and the total po- 
larization Pol = (iV| - iV|)/(7V| -I- iV|) for the trapped 
Fermi gas. Here N„ is obtained by integrating the local 
density n^ir) over the trapped volume. These are di- 
rectly measurable quantities in the experiments (see Fig 
Eland the reference [l^])- In our result, we observe that 
the radial boundary of the first order phase transition 
appears at Pol lower than the critical Pole ~ 0.65. Pole 
is defined when the radius ri ^ (see the caption of Fig 
|4]) . This is qualitatively different from the value given by 
the mean field method 24] where at the unitary regime 
the Pole is at ^ 1 and the superfluid core appears at 
any Pol < 1. Also, the 7 « 0.85 from the same mean 
field work in larger discrepancy with the experiments. 
Our Pole is somewhat lower than those given by the ref- 
erences and m {Pole « 0.77 and 0.75 respectively). 
The differences can be attributed to the calculation meth- 
ods and the experimental conditions. 

In conclusion, we implemented a fully ab initio method 
for calculating the equation of state of the unitary Fermi 
gas at zero temperature. The sign problem of the spin 
imbalanced systems makes the Monte Carlo integration 
somewhat inefficient but still possible for the and r 
considered. The comparisons at the extremes of the = 
Ni and TV^ = 1 with the available literature produce good 
match. This gave us confidence that in < a; < I, our 
result is close to the accurate equation of state. We were 
also able to extract the actual yi as the limit of y{x) 
at 2: — !■ 1. By using the thermodynamic relations and 
LDA we could draw the local densities of the trapped gas 
where we can locate the phase transition radii at different 
total polarizations (Pol). The behavior of the pressure in 
different quantum phases was also studied. This method 
is general and we could easily extend to regimes off the 
unitary. 
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